CLL Cytokine Screen 2021: Figure 5
19 July 2021
In this sub-vignette we present the analysis and source code for figure 5. This sub-vignette can be built along with all other sub-vignettes, by running CLLCytokineScreen2021.Rmd.
1 Set up
Load Libraries
library(patchwork)
library(ggplot2)
library(data.table)
library(magrittr)
library(ggplot2)
library(pheatmap)
library(tidyr)
library(dplyr)Set plot directory
plotDir = ifelse(exists(".standalone"), "", "../../inst/figs/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)2 Load data
Raw
#df: tibble containing all screening viability data
load( "../../data/df.RData")#from tsvs
df <- read.table(file= "../../inst/extdata/df.txt",header = TRUE) %>% as_tibble()
df$Cytokine <- factor(df$Cytokine, levels= c("No Cytokine",
"Resiquimod",
"IL-4",
"TGF-b1",
"IL-1b",
"Interferon gamma",
"SDF-1a",
"sCD40L" ,
"sCD40L+IL-4",
"soluble anti-IgM",
"CpG ODN",
"IL-6",
"IL-10",
"IL-21",
"HS-5 CM",
"IL-15",
"BAFF",
"IL-2" ))
df$Drug <- factor(df$Drug, levels =c("DMSO",
"BAY-11-7085",
"Everolimus",
"Fludarabine",
"I-BET 762",
"Ibrutinib","Idelalisib",
"Luminespib",
"Nutlin-3a",
"PRT062607",
"Pyridone-6",
"Ralimetinib",
"Selumetinib"))Process data
# List of Cytokines, Drugs and treatments
##Drugs
thedrugs <- unique(df$Drug)
##Cytokines
thecytokines <- unique(df$Cytokine)
#Treatments
#drugs
drugtreatments <- list()
drugtreatments <-
lapply(thedrugs, function(x){
paste("treatment_drug", x, sep='')
})
#cytokines
cytokinetreatments <- list()
cytokinetreatments <-
lapply(thecytokines, function(y){
paste("treatment_cytokine", y, sep='')
})
# Tidy up data frame to use in linear model
df <-
dplyr::filter(df,
#use high concentrations only
Drug_Concentration %in% c("High","None")) %>%
#select columns for linear model
dplyr::select(PatientID, Drug, Cytokine, Log) %>%
#rename columns
plyr::rename(c("Drug"="treatment_drug", "Cytokine"="treatment_cytokine", "Log"="Viability"))3 Set aesthetics
source("../../R/themes_colors.R")4 Pre - work
Fit the linear model
# define the base level per treatment_type as the "no"-treatment
df$treatment_drug <- as.factor(df$treatment_drug)
df$treatment_cytokine <- as.factor(df$treatment_cytokine)
df$treatment_drug %<>% relevel("DMSO")
df$treatment_cytokine %<>% relevel("No Cytokine")
# inspect design matrix of interaction model
# model.matrix(~treatment_drug + treatment_cytokine + treatment_drug:treatment_cytokine, df)
# fit model with interaction
fit <- lm(Viability ~ treatment_drug * treatment_cytokine, df)Extract p values
#extract p values
pvalues <- summary(fit)$coefficients[,4]
#order as a dataframe
pvaldf <- as.data.frame(as.matrix(pvalues)) %>%
setDT(keep.rownames = TRUE)
#rename columns
colnames(pvaldf) <- c("treatment", "pvalue")
#filter out single agent treatments
singletreatments <- c(drugtreatments, cytokinetreatments, "(Intercept)")
pvaldf <- dplyr::filter(pvaldf, !treatment %in% singletreatments)
pvaldf <-
pvaldf %>%
#remove "treatment" prefix
mutate_at(vars(treatment), funs(as.character(gsub("treatment_drug", "", .)))) %>%
mutate_at(vars(treatment), funs(as.character(gsub("treatment_cytokine", "", .))))## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE)) #split drug:cytokine into two columns
pvaldf <-
data.frame(pvaldf, do.call(rbind, strsplit(pvaldf$treatment, split = ":", fixed = TRUE)))
#renove treatment column
pvaldf <- pvaldf[,c("pvalue","X1", "X2")]
#rename columns
colnames(pvaldf) <- c("pvalue","drug", "Cytokine")5 Plot Figures
5.1 Fig 5A
Dummy plots
dummydata <-
as.data.frame(list(c(1,2,3,4),
c(0,0,0,0),
c(-0.6,-0.6,0,-0.1),
c(0.6,0.6,0.1,-0.1),
c(0.5,-0.5,0.8,-0.8)),
col.names=c("plot",
"coord1",
"coord2",
"coord3",
"coord4")) %>%
as_tibble()%>%
mutate(expected = coord2 + coord3, measured = coord4) %>%
pivot_longer(cols = coord1:measured, names_to = c("plot"), values_to="y.value", names_repair="unique" ) %>%
dplyr::rename(Plot = plot...1, x.value = plot...2) %>%
mutate(x.value = factor(x.value, levels = c("coord1", "coord2", "coord3", "coord4", "expected", "measured"))) %>%
#add annotation of interaction type
mutate(Plot = factor(Plot,
labels = c("stimulus inhibits drug",
"drug inhibits stimulus",
"synergistic pro-survival effect",
"synergistic toxicity")))## New names:
## * plot -> plot...1
## * plot -> plot...2gg_dummy =
#for each interaction type, make a dummy plot
lapply(unique(dummydata$Plot), function(i){
tmp <- dplyr::filter(dummydata, Plot==i )
#get value for interaction
interaction =
unlist(dplyr::select(dplyr::filter(tmp, x.value=="measured"), y.value)) -
unlist(dplyr::select(dplyr::filter(tmp,x.value=="expected"), y.value))
#plot dummy plot
gg_dummy <-
ggplot(dplyr::filter(tmp, grepl('coord', x.value)),
aes(x=x.value, y=y.value, group=Plot))+
geom_hline(yintercept = 0, linetype=2, color="black")+
geom_point(colour=ifelse(interaction > 0, "#A6093D","#003DA5"), size=1)+
geom_line(colour=ifelse(interaction > 0, "#A6093D","#003DA5"), size=1)+
geom_segment(data=dplyr::filter(tmp, grepl('expected', x.value)), aes(x = 3.5, y = y.value, xend = 4.5, yend = y.value), color="blue", size=2) +
geom_segment(data=dplyr::filter(tmp, grepl('measured', x.value)), aes(x = 3.5, y = y.value, xend = 4.5, yend = y.value), color="black", size=2)+
theme_bw() +
t1 +
xlab("") +
ylab("") +
coord_cartesian(ylim = c(-1, 1)) +
scale_x_discrete(labels=c( "coord1"="DMSO", "coord2"="Drug","coord3"="Stimulus","coord4"="Drug +\nStimulus")) +
theme(axis.ticks.x= element_blank(),axis.ticks.y = element_blank(),
plot.tag=element_text(size = 20)) +
scale_y_continuous(breaks = c(-1,0,1), sec.axis =
sec_axis(trans=~exp(1)^.*100,
breaks=c(50,100,200),
labels=c("50%", "100%", "200%")))
})
#construct plot
Fig5A <- wrap_elements(
gg_dummy[[1]] +
labs(title=expression(paste("Positive ", beta["int"] ))) +
ylab("Antagonism")+
gg_dummy[[2]] +
labs(title=expression(paste("Negative ", beta["int"] )))+
gg_dummy[[3]] +
ylab("Synergism") +
gg_dummy[[4]]+
plot_layout(nrow=2) +
plot_annotation(tag_levels = "I")
)
Fig5A
5.2 Fig 5B
Barplot quantifying interaction types
fit_as_df <-
fit$coefficients %>%
as.matrix() %>%
as.data.frame() %>%
setDT(keep.rownames = TRUE) %>%
as_tibble() %>%
dplyr::rename(Condition = rn, comb_coefficient = V1) %>%
dplyr::filter(grepl("drug", Condition)&grepl("cytokine", Condition)) %>%
mutate(Condition = gsub("treatment_drug", "", Condition)) %>%
mutate(Condition = gsub("treatment_cytokine", "", Condition)) %>%
separate(Condition, c("Drug", "Cytokine"), sep=":")
drug_fit_as_df <-
fit$coefficients %>%
as.matrix() %>%
as.data.frame() %>%
setDT(keep.rownames = TRUE) %>%
as_tibble() %>%
dplyr::rename(Drug = rn, drug_coefficent=V1) %>%
dplyr::filter(grepl("drug", Drug)&!grepl("cytokine", Drug)) %>%
mutate(Drug = gsub("treatment_drug", "", Drug))
cyt_fit_as_df <-
fit$coefficients %>%
as.matrix() %>%
as.data.frame() %>%
setDT(keep.rownames = TRUE) %>%
as_tibble() %>%
dplyr::rename(Cytokine=rn, cyt_coefficent=V1) %>%
dplyr::filter(!grepl("drug", Cytokine)&grepl("cytokine", Cytokine)) %>%
mutate(Cytokine = gsub("treatment_cytokine", "", Cytokine))
fit_combinations <-
left_join(fit_as_df, drug_fit_as_df, by = "Drug") %>%
left_join(cyt_fit_as_df, by = "Cytokine") %>%
mutate(predicted_coefficient=cyt_coefficent+drug_coefficent) %>%
mutate(obs_coef=predicted_coefficient+comb_coefficient) %>%
left_join(mutate(pvaldf,drug=as.character(drug),Cytokine=as.character(Cytokine)), by=c("Drug"="drug", "Cytokine"="Cytokine")) %>%
dplyr::filter(pvalue<=0.05) %>%
mutate(Category_Symbol=case_when(
comb_coefficient>=0&(obs_coef<cyt_coefficent|obs_coef<drug_coefficent)~ "I",
comb_coefficient>=0&(obs_coef>cyt_coefficent&obs_coef>drug_coefficent)~ "III",
comb_coefficient<0&(obs_coef>cyt_coefficent|obs_coef>drug_coefficent)~ "II",
comb_coefficient<0&(obs_coef<cyt_coefficent&obs_coef<drug_coefficent)~ "IV"
)) %>%
mutate(Category=case_when(
Category_Symbol=="I"~ "Positive Antagonistic",
Category_Symbol=="III"~ "Positive Synergistic",
Category_Symbol=="II"~ "Negative Antagonistic",
Category_Symbol=="IV"~ "Negative Synergistic"
)) %>%
mutate( Category=factor(Category, levels=c("Positive Antagonistic","Negative Antagonistic","Positive Synergistic","Negative Synergistic")))
Fig5B <-
fit_combinations %>%
dplyr::group_by(Category) %>%
tally() %>%
ggplot(aes(x=Category, y=n, fill=Category))+
geom_bar(width = 1, stat = "identity")+
t1+
guides(fill="none")+
xlab("")+ ylab("")+
scale_fill_manual(values = c("#A6093D","#003DA5","#A6093D","#003DA5"))
Fig5B
5.3 Fig 5C
Heatmap of significant drug : stimulus interaction coefficients
Get interaction coefficients matrix for plotting heatmap
#extract beta interaction values
betavalues <- summary(fit)$coefficients[,1]
#create matrix of beta values
betadf <- as.data.frame(as.matrix(betavalues)) %>%
setDT(keep.rownames = TRUE)
#rename columns
colnames(betadf) <- c("treatment", "betavalue")
#remove beta values for drugs and cytokines alone
betadf <- dplyr::filter(betadf, !treatment %in% singletreatments)
#remove treatment prefix
betadf <- betadf %>%
mutate_at(vars(treatment), funs(as.character(gsub("treatment_drug", "", .)))) %>%
mutate_at(vars(treatment), funs(as.character(gsub("treatment_cytokine", "", .))))
#split drug:cytokine by colon
betadf <- data.frame(betadf, do.call(rbind, strsplit(betadf$treatment, split = ":", fixed = TRUE)))
#select columns of interest
betadf <- betadf[,c("betavalue","X1", "X2")]
#rename columns
colnames(betadf) <- c("betavalue","drug", "Cytokine")#set significance
a <- 0.05
#make a matrix of beta values, where beta = 0 if corresponding p val is > than significance threshold
pvalmat <- xtabs(pvalue~drug+Cytokine, data=pvaldf)
bvalmat <- xtabs(betavalue~drug+Cytokine, data=betadf)
#check in same order
bvalmat[,colnames(pvalmat)] ## Cytokine
## drug BAFF CpG ODN HS-5 CM IL-10
## BAY-11-7085 -5.604631e-02 -3.181140e-01 1.492971e-03 -1.575975e-02
## Everolimus 3.655902e-02 1.497761e-02 5.377161e-02 7.948169e-02
## Fludarabine -2.337342e-02 1.903315e-01 6.711677e-02 -3.913832e-02
## I-BET 762 4.216039e-02 -2.331726e-01 -7.116831e-02 -2.652204e-02
## Ibrutinib -2.625863e-02 -2.269246e-01 7.477386e-02 2.731021e-02
## Idelalisib 8.275484e-03 -2.652441e-01 3.386327e-02 3.058002e-02
## Luminespib -2.482555e-01 -2.591311e-01 -3.940346e-01 -6.458359e-02
## Nutlin-3a -1.258362e-03 1.577514e-01 7.700330e-02 9.299536e-02
## PRT062607 -6.653378e-02 -8.644230e-02 1.557073e-02 -8.013149e-03
## Pyridone-6 -8.023060e-04 -8.100430e-03 2.925389e-02 3.243689e-02
## Ralimetinib -3.093962e-02 -1.717127e-01 1.858811e-02 -1.073542e-02
## Selumetinib 5.947874e-03 -1.304805e-01 7.060345e-02 2.145764e-02
## Cytokine
## drug IL-15 IL-1b IL-2 IL-21
## BAY-11-7085 -5.816512e-03 2.152596e-02 -2.134484e-02 2.072687e-03
## Everolimus 2.725033e-02 3.500697e-02 -2.677401e-03 3.495098e-02
## Fludarabine -4.310581e-02 3.505833e-02 -2.214672e-02 -4.860353e-02
## I-BET 762 9.141791e-04 2.740045e-03 -2.907140e-02 -2.523721e-02
## Ibrutinib 5.752202e-02 1.849707e-03 4.678019e-02 3.710308e-02
## Idelalisib 5.658295e-02 1.936971e-02 4.584077e-02 3.477902e-02
## Luminespib -6.061613e-02 -1.263584e-01 -6.993392e-02 -7.742445e-02
## Nutlin-3a 1.046024e-01 2.482328e-02 1.463441e-01 7.912615e-02
## PRT062607 -3.957221e-03 -2.216149e-02 5.916628e-03 8.482059e-03
## Pyridone-6 5.308392e-02 1.514414e-02 3.040022e-02 -7.237789e-04
## Ralimetinib 1.202065e-02 -2.568524e-03 -1.468196e-02 -4.467606e-02
## Selumetinib 5.942633e-02 1.780087e-02 3.092410e-02 2.662341e-02
## Cytokine
## drug IL-4 IL-6 Interferon gamma Resiquimod
## BAY-11-7085 2.235761e-02 4.554573e-02 5.434015e-02 -3.639915e-01
## Everolimus 1.332949e-01 6.817609e-03 5.469766e-02 2.765104e-02
## Fludarabine 3.953363e-01 -2.469675e-02 4.886592e-02 2.265297e-01
## I-BET 762 -2.665207e-02 -2.262955e-02 1.357931e-02 -3.265356e-01
## Ibrutinib 1.836975e-01 6.642272e-02 1.101352e-01 -2.509692e-01
## Idelalisib 1.657969e-01 6.363235e-02 1.171099e-01 -2.064261e-01
## Luminespib -4.429811e-02 -3.220023e-03 -2.221338e-02 -7.116024e-02
## Nutlin-3a 5.337703e-01 1.136707e-01 1.433726e-01 1.426722e-01
## PRT062607 1.388865e-01 -1.097500e-02 2.337567e-02 -3.612143e-01
## Pyridone-6 -2.305619e-02 6.847247e-02 6.271102e-02 1.349796e-03
## Ralimetinib 7.075108e-02 2.282944e-03 2.190150e-01 -1.118603e-01
## Selumetinib 1.174186e-01 6.322721e-02 8.349273e-02 -1.148979e-01
## Cytokine
## drug sCD40L sCD40L+IL-4 SDF-1a soluble anti-IgM
## BAY-11-7085 2.049505e-02 -1.619583e-01 2.831627e-02 -6.597054e-03
## Everolimus 5.178276e-02 1.386887e-01 -2.677009e-03 5.541420e-04
## Fludarabine 1.020548e-02 3.687631e-01 5.075778e-02 3.435734e-02
## I-BET 762 4.068479e-02 6.576262e-03 -4.944505e-03 -2.502231e-02
## Ibrutinib -4.740237e-03 8.873182e-02 2.038604e-02 1.260542e-02
## Idelalisib 2.027841e-02 8.430112e-02 -2.252409e-05 2.338468e-02
## Luminespib -1.543305e-01 -3.482382e-01 -8.004006e-02 -4.275447e-01
## Nutlin-3a 1.834396e-02 4.643749e-01 1.258999e-01 5.912671e-03
## PRT062607 -8.838768e-02 1.457566e-01 -2.877775e-02 -6.326385e-02
## Pyridone-6 3.934846e-02 -1.855970e-01 3.145554e-02 -6.082148e-03
## Ralimetinib -1.691833e-04 -3.878371e-02 1.575747e-03 -2.273022e-02
## Selumetinib 5.122307e-02 4.874905e-02 -2.510831e-02 2.421239e-02
## Cytokine
## drug TGF-b1
## BAY-11-7085 -4.509014e-03
## Everolimus -8.310092e-03
## Fludarabine -4.636054e-02
## I-BET 762 -9.396045e-02
## Ibrutinib 1.958423e-02
## Idelalisib 9.050594e-03
## Luminespib -8.478118e-02
## Nutlin-3a -1.529464e-02
## PRT062607 2.369375e-03
## Pyridone-6 3.081690e-02
## Ralimetinib 1.357279e-02
## Selumetinib 2.234705e-02bvalmat[rownames(pvalmat),]## Cytokine
## drug BAFF CpG ODN HS-5 CM IL-10
## BAY-11-7085 -5.604631e-02 -3.181140e-01 1.492971e-03 -1.575975e-02
## Everolimus 3.655902e-02 1.497761e-02 5.377161e-02 7.948169e-02
## Fludarabine -2.337342e-02 1.903315e-01 6.711677e-02 -3.913832e-02
## I-BET 762 4.216039e-02 -2.331726e-01 -7.116831e-02 -2.652204e-02
## Ibrutinib -2.625863e-02 -2.269246e-01 7.477386e-02 2.731021e-02
## Idelalisib 8.275484e-03 -2.652441e-01 3.386327e-02 3.058002e-02
## Luminespib -2.482555e-01 -2.591311e-01 -3.940346e-01 -6.458359e-02
## Nutlin-3a -1.258362e-03 1.577514e-01 7.700330e-02 9.299536e-02
## PRT062607 -6.653378e-02 -8.644230e-02 1.557073e-02 -8.013149e-03
## Pyridone-6 -8.023060e-04 -8.100430e-03 2.925389e-02 3.243689e-02
## Ralimetinib -3.093962e-02 -1.717127e-01 1.858811e-02 -1.073542e-02
## Selumetinib 5.947874e-03 -1.304805e-01 7.060345e-02 2.145764e-02
## Cytokine
## drug IL-15 IL-1b IL-2 IL-21
## BAY-11-7085 -5.816512e-03 2.152596e-02 -2.134484e-02 2.072687e-03
## Everolimus 2.725033e-02 3.500697e-02 -2.677401e-03 3.495098e-02
## Fludarabine -4.310581e-02 3.505833e-02 -2.214672e-02 -4.860353e-02
## I-BET 762 9.141791e-04 2.740045e-03 -2.907140e-02 -2.523721e-02
## Ibrutinib 5.752202e-02 1.849707e-03 4.678019e-02 3.710308e-02
## Idelalisib 5.658295e-02 1.936971e-02 4.584077e-02 3.477902e-02
## Luminespib -6.061613e-02 -1.263584e-01 -6.993392e-02 -7.742445e-02
## Nutlin-3a 1.046024e-01 2.482328e-02 1.463441e-01 7.912615e-02
## PRT062607 -3.957221e-03 -2.216149e-02 5.916628e-03 8.482059e-03
## Pyridone-6 5.308392e-02 1.514414e-02 3.040022e-02 -7.237789e-04
## Ralimetinib 1.202065e-02 -2.568524e-03 -1.468196e-02 -4.467606e-02
## Selumetinib 5.942633e-02 1.780087e-02 3.092410e-02 2.662341e-02
## Cytokine
## drug IL-4 IL-6 Interferon gamma Resiquimod
## BAY-11-7085 2.235761e-02 4.554573e-02 5.434015e-02 -3.639915e-01
## Everolimus 1.332949e-01 6.817609e-03 5.469766e-02 2.765104e-02
## Fludarabine 3.953363e-01 -2.469675e-02 4.886592e-02 2.265297e-01
## I-BET 762 -2.665207e-02 -2.262955e-02 1.357931e-02 -3.265356e-01
## Ibrutinib 1.836975e-01 6.642272e-02 1.101352e-01 -2.509692e-01
## Idelalisib 1.657969e-01 6.363235e-02 1.171099e-01 -2.064261e-01
## Luminespib -4.429811e-02 -3.220023e-03 -2.221338e-02 -7.116024e-02
## Nutlin-3a 5.337703e-01 1.136707e-01 1.433726e-01 1.426722e-01
## PRT062607 1.388865e-01 -1.097500e-02 2.337567e-02 -3.612143e-01
## Pyridone-6 -2.305619e-02 6.847247e-02 6.271102e-02 1.349796e-03
## Ralimetinib 7.075108e-02 2.282944e-03 2.190150e-01 -1.118603e-01
## Selumetinib 1.174186e-01 6.322721e-02 8.349273e-02 -1.148979e-01
## Cytokine
## drug sCD40L sCD40L+IL-4 SDF-1a soluble anti-IgM
## BAY-11-7085 2.049505e-02 -1.619583e-01 2.831627e-02 -6.597054e-03
## Everolimus 5.178276e-02 1.386887e-01 -2.677009e-03 5.541420e-04
## Fludarabine 1.020548e-02 3.687631e-01 5.075778e-02 3.435734e-02
## I-BET 762 4.068479e-02 6.576262e-03 -4.944505e-03 -2.502231e-02
## Ibrutinib -4.740237e-03 8.873182e-02 2.038604e-02 1.260542e-02
## Idelalisib 2.027841e-02 8.430112e-02 -2.252409e-05 2.338468e-02
## Luminespib -1.543305e-01 -3.482382e-01 -8.004006e-02 -4.275447e-01
## Nutlin-3a 1.834396e-02 4.643749e-01 1.258999e-01 5.912671e-03
## PRT062607 -8.838768e-02 1.457566e-01 -2.877775e-02 -6.326385e-02
## Pyridone-6 3.934846e-02 -1.855970e-01 3.145554e-02 -6.082148e-03
## Ralimetinib -1.691833e-04 -3.878371e-02 1.575747e-03 -2.273022e-02
## Selumetinib 5.122307e-02 4.874905e-02 -2.510831e-02 2.421239e-02
## Cytokine
## drug TGF-b1
## BAY-11-7085 -4.509014e-03
## Everolimus -8.310092e-03
## Fludarabine -4.636054e-02
## I-BET 762 -9.396045e-02
## Ibrutinib 1.958423e-02
## Idelalisib 9.050594e-03
## Luminespib -8.478118e-02
## Nutlin-3a -1.529464e-02
## PRT062607 2.369375e-03
## Pyridone-6 3.081690e-02
## Ralimetinib 1.357279e-02
## Selumetinib 2.234705e-02bvalmat[pvalmat >= a] = 0
#transform matrix
bvalmat <- t(bvalmat)
tree <-
pheatmap(bvalmat, silent = TRUE)
cyt_order <- rownames(bvalmat[tree$tree_row[["order"]],])
drug_order <- colnames(bvalmat[,tree$tree_col[["order"]]])
Fig5C <-
#append p values
left_join(betadf, pvaldf, by = c("drug", "Cytokine")) %>%
#make beta value 0 if p value is < sig
mutate(betavalue = ifelse(pvalue<a, betavalue, 0)) %>%
#put drugs in order of clustered heatmap
mutate(drug=factor(drug, levels = drug_order)) %>%
#put stimuli in order of clustered heatmap
mutate(Cytokine=factor(Cytokine, levels = rev(cyt_order))) %>%
#make heatmap with ggplot
ggplot(aes(x=drug, y=Cytokine))+
geom_tile(aes(fill=betavalue),color = "grey")+
#set colour scale
scale_fill_gradientn(colors=c("#003DA5", "white", "#A6093D"), limits=c(-0.6,0.6))+
#set theme
t2+
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank(),
axis.ticks.x =element_blank(),
panel.background = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
legend.key.height=unit(2.5, "cm"))+
labs(fill = expression(paste(beta["int"], "-value")))+
#add category
geom_text(data = fit_combinations, aes(x=Drug, label=Category_Symbol) )+
scale_y_discrete(labels=c("TGF-b1"="TGF-\u03B21", "sCD40L+IL-4"="sCD40L + IL4", "IL-1b"="IL1\u03B2", "IL-4"="IL4", "IL-6"="IL6","IL-15"="IL15","IL-10"="IL10", "IL-21"="IL21","IL-2"="IL2", "Interferon gamma"= "Interferon \u03B3", "SDF-1a"="SDF-1\u03B1"))
Fig5C
5.4 Fig 5D - I
Preparation: Line plots of examples of drug and stimulus combinations that show an interaction
#get lists of drugs and cytokines without baseline treatments
thedrugs.only <- thedrugs %>% setdiff("DMSO") %>% sort()
thecytokines.only <- thecytokines %>% setdiff("No Cytokine")%>% sort()
# get a list of interactions for which p value of beta interaction is significant
sign_conditions <-
pvaldf %>%
dplyr::filter(pvalue<=0.05) %>%
mutate(condition = paste0("treatment_drug:",
drug,
":treatment_cytokine:",
Cytokine))
gg <- vector(mode="list", length=length(sign_conditions$condition))
names(gg) <- sign_conditions$condition
plotList <-
lapply(names(gg), function(i){
#preparations for the interaction plot
drug <- dplyr::filter(sign_conditions, condition ==i)$drug
drugTreat <- paste("treatment_drug", drug, sep='')
cyt <- dplyr::filter(sign_conditions, condition ==i)$Cytokine
cytokineTreat <-paste("treatment_cytokine", cyt, sep='')
Category <- fit_combinations[which(fit_combinations$Drug==drug &
fit_combinations$Cytokine==cyt),]$Category
#extract coefficients to plot
## baseline effect
baseline <- fit$coefficients["(Intercept)"]
## single treatment effect for type_1 = "a"
single_1 <- fit$coefficients[drugTreat]
## single treatment effect for type_2 = "b"
single_2 <- fit$coefficients[cytokineTreat]
# interaction Viability
interaction <- fit$coefficients[paste(drugTreat, cytokineTreat, sep=':')]
## observed single treatment values
observed_single_1 <- baseline + single_1
observed_single_2 <- baseline + single_2
## observed double treatment value
observed_double <- baseline + single_1 + single_2 + interaction
## predicted double treatment value (based only on single effects)
predicted_double <- baseline + single_1 + single_2
#prepare labels and aesthetics for plotting
# y_label <- case_when(cyt=="Interferon gamma"~"Interferon\ngamma",
# cyt=="soluble anti-IgM"~"soluble\nanti-IgM",
# TRUE~cyt)
xlabels <- c("DMSO", drug,
case_when(cyt=="IL-4"~"IL4",
cyt=="Interferon gamma"~"Interferon \u03B3",
cyt=="sCD40L+IL-4"~"sCD40L + IL4",
cyt=="soluble anti-IgM"~"soluble\nanti-IgM",
TRUE~cyt),
paste(drug, " +\n",
case_when(cyt=="IL-4"~"IL4",
cyt=="Interferon gamma"~"Interferon \u03B3",
cyt=="sCD40L+IL-4"~"sCD40L + IL4",
TRUE~cyt),
sep=""))
segment_size <- 2
plotTab <-
df %>%
dplyr::filter(treatment_drug %in% c("DMSO" , drug)) %>%
dplyr::filter(treatment_cytokine %in% c("No Cytokine", cyt)) %>%
dplyr::mutate(treatment_combination =
#if baseline treatment
ifelse(treatment_drug=="DMSO" &
treatment_cytokine=="No Cytokine",
"DMSO",
#if drug only treatment
ifelse(treatment_drug == drug &
treatment_cytokine =="No Cytokine",
paste0("Drug=", drug),
#if cytokine only treatment
ifelse(treatment_drug == "DMSO" &
treatment_cytokine == cyt,
paste0("Cytokine=", cyt),
#otherwise combinatorial treatment
paste0("Drug=", drug, "\nCytokine=", cyt))))) %>%
#make treatment_combination as a factor
mutate(treatment_combination = factor(treatment_combination,
levels=c("DMSO",
paste0("Drug=", drug),
paste0("Cytokine=", cyt),
paste0("Drug=", drug, "\nCytokine=", cyt)
)))
#make plot
ggplot(plotTab, aes(treatment_combination, Viability, group=PatientID)) +
geom_hline(yintercept = 0, linetype=2, color="black")+
lemon::geom_pointline( size=1,
colour=ifelse(interaction > 0,
"#A6093D","#003DA5"), alpha=0.3) +
#add predicted viability with control treatment
geom_segment( size=segment_size,
aes(x=1-.2, xend=1+.2, y=baseline, yend=baseline),
colour = "black") +
#add predicted viability with drug only treatment
geom_segment( size=segment_size,
aes(x=2-.2, xend=2+.2, y=observed_single_1,yend=observed_single_1),
colour="black") +
#add predicted viability with cytokine only treatment
geom_segment( size=segment_size,
aes(x=3-.2, xend=3+.2, y=observed_single_2, yend=observed_single_2),
colour="black") +
#add predicted viability with both treatments
geom_segment( size=segment_size,
aes(x=4-.2, xend=4+.2, y=observed_double, yend=observed_double),
colour="black") +
#add predicted viability with both treatments, without interaction
geom_segment( size=segment_size,
aes(x=4-.2, xend=4+.2, y=predicted_double,
yend=predicted_double), color="blue") +
#add category and combination as title
ggtitle(paste( Category,": ", drug, " & ",
case_when(cyt=="IL-4"~"IL4",
cyt=="Interferon gamma"~"Interferon \u03B3",
cyt=="sCD40L+IL-4"~"sCD40L + IL4",
TRUE~cyt), sep="")) +
#add theme
t2 + theme(axis.text.x = element_text(size=fontsize+5)) +
theme(plot.tag=element_text(size = 30))+
scale_x_discrete(labels = xlabels) +
ylab("Log (Viability)") +
xlab("")
} )
names(plotList) <- names(gg)5.5 Fig 5D
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:IL-4`)
5.6 Fig 5E
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:Interferon gamma`)
5.7 Fig 5F
wrap_plots(plotList$`treatment_drug:Pyridone-6:treatment_cytokine:sCD40L+IL-4`)
5.8 Fig 5G
wrap_plots(plotList$`treatment_drug:Ralimetinib:treatment_cytokine:Interferon gamma`)
5.9 Fig 5H
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:CpG ODN`)
5.10 Fig 5G
wrap_plots(plotList$`treatment_drug:Luminespib:treatment_cytokine:soluble anti-IgM`)
6 Arrange plots
design1 <-"
1223
4455
6677
8899
"
Fig5A + theme(plot.tag=element_text(size = 30)) +
wrap_elements(Fig5B) + theme(plot.tag=element_text(size = 30)) +
wrap_elements(Fig5C) + theme(plot.tag=element_text(size = 30)) +
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:IL-4`)+
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:Interferon gamma`)+
wrap_plots(plotList$`treatment_drug:Pyridone-6:treatment_cytokine:sCD40L+IL-4`)+
wrap_plots(plotList$`treatment_drug:Ralimetinib:treatment_cytokine:Interferon gamma`)+
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:CpG ODN`)+
wrap_plots(plotList$`treatment_drug:Luminespib:treatment_cytokine:soluble anti-IgM`)+
plot_layout(design=design1,
heights = c(1,0.5,0.5,0.5),
widths = c(0.92,0.08,0.05,0.95)) +
plot_annotation( tag_levels = "A", theme=theme(plot.title = element_text(face="plain", size=28, hjust = 0.5)))
7 Appendix
Sys.info()## sysname
## "Darwin"
## release
## "19.6.0"
## version
## "Darwin Kernel Version 19.6.0: Thu May 6 00:48:39 PDT 2021; root:xnu-6153.141.33~1/RELEASE_X86_64"
## nodename
## "mac-huber20.academic.christs.cam.ac.uk.pool.embl.de"
## machine
## "x86_64"
## login
## "holly"
## user
## "holly"
## effective_user
## "holly"sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.7 tidyr_1.1.3 pheatmap_1.0.12 magrittr_2.0.1
## [5] data.table_1.14.0 ggplot2_3.3.5 patchwork_1.1.1 BiocStyle_2.20.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.6 highr_0.9 plyr_1.8.6
## [4] RColorBrewer_1.1-2 pillar_1.6.1 compiler_4.1.0
## [7] BiocManager_1.30.16 tools_4.1.0 digest_0.6.27
## [10] lattice_0.20-44 evaluate_0.14 lifecycle_1.0.0
## [13] tibble_3.1.2 gtable_0.3.0 pkgconfig_2.0.3
## [16] rlang_0.4.11 DBI_1.1.1 magick_2.7.2
## [19] yaml_2.2.1 lemon_0.4.5 xfun_0.24
## [22] gridExtra_2.3 withr_2.4.2 stringr_1.4.0
## [25] knitr_1.33 generics_0.1.0 vctrs_0.3.8
## [28] grid_4.1.0 tidyselect_1.1.1 glue_1.4.2
## [31] R6_2.5.0 fansi_0.5.0 rmarkdown_2.9
## [34] bookdown_0.22.3 farver_2.1.0 purrr_0.3.4
## [37] scales_1.1.1 ellipsis_0.3.2 htmltools_0.5.1.1
## [40] assertthat_0.2.1 colorspace_2.0-2 labeling_0.4.2
## [43] utf8_1.2.1 stringi_1.6.2 munsell_0.5.0
## [46] crayon_1.4.1